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BACKGROUND OF THE INVENTION 



Iff"" Technical Field 

W 

111 The invention relates to the field of molecular similarity searching, and, more 

*S specifically, similarity searching in databases of three dimensional molecular structures. 

Description of the Related Art 

In the field of drug design, where one is attempting to expand the number of lead 
1 5 compounds that show activity toward a particular therapeutic target, structural information about 
the target is often lacking or unavailable. Similarity searching in files of chemical compounds is 
a common way to uncover new leads in such situations. Typically, one or more compounds that 
are known to be active toward the target of interest are selected, and a feature scheme is defined 
that characterizes the molecular properties of interest. Features are derived from the selected 
20 structures and used to search against a database of structures that have been keyed under the 
same feature scheme. 
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Feature schemes may be structural (three-dimensional) in nature or topological (derived 
solely from the molecular graph). Features that are three-dimensional (3D) characterize a whole 
or part of a particular conformation of a molecule, and thus are dependent on the particular 
conformations of the molecules stored in the database. 3D features can include a) 
pharmacophore descriptors, such as distances, angles, or dihedral angular relationships between 
key groups (see Martin, Y. et. al., A fast new approach to pharmacophore mapping and its 
application to dopaminergic and benzodiazepine agonists. J. Comput. -Aided Mol. Des. 1992, Vol 
6, pp. 475-486; b) Surfaces characterizations (see Perkins et. al. Molecular surface-volume and 
property matching to superpose flexible dissimilar molecules. J. Comput.-Aided Mol. Des. 1995, 
Vol 9, 479-490) ; or c) Field-based properties that characterize regions of a molecule (see Willet. 
et. al. in Similarity searching in files of three-dimensional chemical structures: Flexible 
field-based searching of molecular electrostatic potentials. J. Chem. Inf. Comput. Sci. 1996, Vol 
36, pp. 900-908). 

Similarity searching for compounds in 3D databases is an important part of lead 
generation, and is commonly practiced in the drug design process (see Klebe G, Structural 
Alignment of Molecules, in 3D QSAR in Drug Design. Theory, Methods, and Practace, and 
Kearsley, S. K. et. al, An alternative method for the alignment of molecular structures: 
Maximizing electrostatic and steric overlap. Tetrahedron Comput. Methodol. 1990, Vol 3, pp. 
615-633). It is useful in expanding the list of active compounds for a therapeutic target, finding 
new uses for existing compounds, getting around a competitors patent, or gaining more insight 
into the nature of the therapeutic target under investigation. There are, however, several ways one 
may conduct such searches, with no one method proven superior or universally applicable. Novel 
procedures are thus a current research interest. 

A common problem that arises in similarity searching is that of preparing an appropriate 
distance metric. The problem arises when one must decide how to weight the relative importance 
of descriptors when evaluating whether two features are similar. The problem is compounded by 
the fact that different contexts warrant different scalings of descriptors. Appropriate distance 
metrics in one context may not be suitable for another. 
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The present invention embodies a novel procedure for a 3D similarity searching that is 
based on the alignment of heuristic property fields. The particular novelty offered by the present 
invention is the independence of the particular property field used, and a context dependent 
scaling procedure that allows a training set to scale the descriptors. 

SUMMARY OF THE INVENTION 

The problems stated above and the related problems of the prior art are addressed with 
the principles of the present invention, similarity searching of molecules based upon statistical 
analysis of descriptor vectors characterizing molecular regions. 

In a training phase, an association criterion is generated by which query regions of a 
query molecule are associated with regions of molecules stored in a database. Preferably, the 
association criterion is based upon statistical analysis of groups of descriptor vectors that 
characterize properties of the regions of the molecules stores in the database. 

In an acquisition phase, for each molecule in a series of molecules, the following steps 
are performed for a given molecule. Data that represents the structure of the given molecule is 
read from persistent memory and used to define a set of three-dimensional regions of space in the 
given molecule. For each region, one or properties of the given molecule are mapped to 
property values for grid points of the region. A multi-map entry is generated that identifies the 
region, and position and orientation of set of axes derived from the property values of the grid 
points of the region. The association criterion generated in the training phase is used generate a 
key for the region, and the entry is stored in the multi-map at a location associated with the key. 

In the recognition phase, data that represents the structure of a query molecule is used to 
define a set of regions in the query molecule. For each region, one or properties of the query 
molecule are mapped to property values for grid points of the query region. The association 
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criterion generated in the training phase is used generate a key for the query region. The 
multi-map entry identified by the key is retrieved and the data stored therein are read from the 
table. For each stored region identified by the retrieved table entry, an hypothesized match is 
constructed and added to a vote table. After processing all of the stored regions identified by the 
retrieved multi-map entry for the set of query regions in the query molecule, one or more entries 
of the vote table is selected, the alignment transformations stored in the selected entries are 
applied to corresponding molecules stored in the database, and the resultant alignment(s) of the 
stored molecule in the query frame is reported to the user via an I/O device. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIGS. 1(A) and 1(B) are block diagrams of computer processing systems wherein the 
methods of the present invention may be embodied. 

FIGS. 2(A) and (B) is a flow chart illustrating the method of the present invention in 
mapping a descriptor vector for an item to a space that optimally discriminates between groups of 
items in accordance with the present invention; 

FIG. 3 is a flow chart illustrating operations of step 205 of FIG. 2 in generating a set of 
component vectors that maximize an F distributed criterion function in accordance with the 
present invention. 

FIG. 4 is a table illustrating a multi-factor design with multi-way analysis of variance. 

FIG. 5 is a flow chart illustrating the training phase of a system that identifies molecules 
within a database of molecules that has similar structure to a query molecule in accordance with 
the present invention. 
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FIGS. 7(A) and (B) is a pictorial illustration of a multi-map, which is part of the system 
that identifies molecules within a database of molecules that has similar structure to a query 
molecule. 

FIGS. 8(A) and (B) is a flow chart illustrating the recognition phase of a system identifies 
molecules within a database of molecules that has similar structure to a query molecule in 
accordance with the present invention. 

FIG. 9 is a flow chart illustrating an exemplary embodiment of the operations in 
constructing the match hypothesis in the recognition phase of the present invention. 

DETAILED DESCRIPTION OF THE INVENTION 

The present invention may be implemented on any computer processing system 
including, for example, a personal computer or a workstation. As shown in FIG. 1(A), a 
computer processing system 100 as may be utilized by the present invention generally comprises 
memory 101, at least one central processing unit (CPU) 103 (one shown), and at least one input 
device 107 (such as a keyboard, mouse, joystick, voice recognition system, or handwriting 
recognition system). In addition, the computer processing system includes a nonvolatile storage 
device 108, such as a ROM or fixed disk drive, that stores an operating system and one or more 
application programs that are loaded into the memory 101 .and executed by the CPU 103. In the 
execution of the operating system and application program(s), the CPU may use data stored in 
the nonvolatile storage device 108 and/or memory 101. In addition, the computer processing 
system includes a graphics adapter 104 coupled between the CPU 103 and a display device 105 
such as a CRT display or LCD display. In addition, the computer processing system may include 
a communication link 109 (such as a network adapter, RF link, or modem) coupled to the CPU 
103 that allows the CPU 103 to communicate with other computer processing systems over the 
communication link, for example over the Internet. The CPU 103 may receive portions of the 
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operating system, portions of the application program(s), or portions of the data used by the CPU 
103 in executing the operating system and application program(s). 

It should be noted that the application program(s) executed by the CPU 103 may perform 
the methods of the present invention described below. Alternatively, portions or all of the 
methods described below may be embodied in hardware that works in conjunction with the 
application program executed by the CPU 103. 

According to the present invention, a computer implemented method provides a mapping 
(i.e., a transformation) for the components of the descriptor vectors for a series of items to a 
space that optimally discriminates between groups of items. With reference to FIGS. 2(A) and 
(B), the operation begins in step 201 wherein a series of items are classified into N groups ( N is 
an integer greater than 1), wherein each group is identified by an identifier / ranging from 1 to 
N. For the sake of description, consider groups that contain m items and descriptor vectors; 
however, the present invention is not limited in this respect and can be used for groups that 
contain a non-uniform number of items and descriptor vectors. In addition, my identifies item j 
belonging to group /'where j ranges from 1 to n h and is the descriptor vector 
corresponding to the item my . Preferably, the descriptor vectors are stored as part of a file in 
persistent storage and loaded into non-persistent storage for use by the CPU 103 as needed. 

In step 203, first data representing co variance between the groups of the items, denoted 
8b , are generated. In addition, second data representing covariance within the items belonging 
to the groups, denoted e w , are generated. Note that both the first data and second data follow a 
chi-square distribution. An example of the operations in determining Sb and e w is provided 

below. In this example, the first data ( 8 b ) have a chi-square distribution with N -1 degrees of 
freedom (where N represents the number of groups of items); and the second data (s w ) have a 

chi-square distribution with 1 - N degrees of freedom (where N represents the number of 
groups of items, rti represents the number of items in group z, and X n t represents the sum of 
fix for the N groups). 
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An example of the operations in determining 8 b and s w is now provided. One may 
ny given d< 
group as follows: 



break any given descriptor vector xy into variation between groups and variations within a 



Xij = a + at + ay , where 

~a represents the mean of all the items (of all the groups); 



a i represents deviation from mean of all the items (a) to the mean of group /; and 

^ -> -» 

U U represents deviation from the mean of group / ( a i ) to the descriptor vector Xy . 



For each group i of items (i.e., i ranging from 1 to N), the sample mean, denoted Xi, , of the 
descriptor vectors within group / is calculated. For example, the sample mean of the descriptor 
vectors in a given group i may be calculated as follows: 

= 7t7 Z where the £ operator sums over the range of j = 1 to 



The sample mean ~a of all the items is calculated. For example, the sample mean ~a may be 
calculated as follows: 



'a -jjT.Xt, where the £ operator sums over the range of i = 1 to N 



For each group i of items (i.e., i ranging from 1 to AO, the deviation at is calculated. For 
example, the deviation at for a given group / may be calculated as follows: 



T> -» 
a, -Xi - a 



For each descriptor vector X// , the deviation fly is calculated. For example, the deviation djj 
may be calculated for a given descriptor vector Xy as follows: 
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The covariance between groups, denoted at , may be calculated as follows: 

£b = X w / <^ ttf , where the 2 operator sums over the range 
of i = 1 to N 

And the covariance within groups, denoted e w , may be calculated as follows: 

&ij •> where the £ X operator sums over the range j = 1 
to ti\ for each group / in the range 1 to N 

Note that the covariance between groups Eb has a chi-square distribution with N -1 degrees of 
freedom (where N represents the number of groups of items); and the covariance within groups 

e w has a chi-square distribution with X ^/ - N degrees of freedom (where N represents the 
number of groups of items, n t represents the number of items in group i, and X #/ represents 
the sum of n t for the N groups). The total covariance, denoted s , is represented as follows: 

8 = £j £j (Xij — a ){Xy — a ) = 6t + £w , where the £ X operator sums over 

the range j = 1 to w,- for each 
group i in the range 1 to N 

In step 205, an F distributed criterion function is used to determine a set of component 
vectors that maximize the criterion function. The criterion function has a numerator and a 
denominator, whereby the numerator is based upon the first data ( 8 b ) generated in step 203 and 
the denominator is based upon the second data ( e w ) generated in step 203. The criterion 
function has the general form: 
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f(w)=c(^4) 



e w i 

where w is a vector in the direction whose variation is 
being tested, and C is a constant based upon the degrees of 
freedom in 8t and e w . 

For example, the constant C may be determined as follows; 

_ 1 /degrees of freedom in £^ \ / (N - \) 

^ ~~ 1 /degrees of freedom in s w 1 / 

where N represents the number of groups of items, 

fit represents the number of items in a group, and X n i 
represents the sum of n, for the N groups, and 

The set of component vectors is determined by solving for those w that maximize the criterion 
function 

f(w) 

. An example of the operation of the computer processing system 100 in 
determining the set of component vectors is set forth below with respect to FIG. 3. For the sake 
of description, the number of component vectors belonging to this set is denoted D. 

In step 207, a loop is performed over the set of component vectors generated in step 205 

(i.e., k = 1 D) to calculate the value, denoted /* , of the criterion function at the given 
component vector w. If the operation set forth below with respect to FIG. 3 is used, the loop of 
step 207 may be performed over the set of vectors "5?, and the value fk for a given vector ~y 
may be calculated as follows: 

/(y)= f (e w )~ J s b (s w )~ J y 



where ~y = (e w ) 2 w , and 

where y is the unit vector corresponding to the vector ~y 
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In addition, in step 207, the value fu of the criterion function /(w) at the given component 
vector w and the associated component vector w (or the value /* of the criterion function f(y) 
at the given vector y and associated vector y) are stored by the processing system 100. 

In step 209, an F distributed statistic is generated for subsets of the set of component 
vectors generated in step 205. The F distributed statistic for a given subset of component vectors 
preferably represents a ratio of the variance between groups of items to the variance within 
groups of items along the given subset of component vectors. In this case, the statistic, denoted 
y/ s , characterizing a given subset of component vectors preferably has following form: 

I, fi 

where fk represents the value of the criterion function 

at a given component vector in the subset, 
C is a constant, 

Ls represents the number of fk values in the given subset 
of component vectors, and 

the 2 operation sums over the Ls fk values in the given 
subset of component vectors. 

Note that in the example above where C = ^ ~~ the V s 's are F distributed with 

r 1/(2, nj-N) 

the ((N-\)Ls) degrees of freedom in the numerator and - N) degrees of freedom in 

the denominator. 



An optimal subset S of the component vectors is then selected by identifying a subset of 
component vectors such that a probability value for the statistic y/ s associated with the subset 
(preferably, the probability value represents that the probability that the statistic y/ s for the 
subset could have been larger by chance) satisfies a predetermined criterion (significance level). 
This significance level represents a threshold at which an hypothesis that the aggregate 
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F-distributed ratios (for the subset) representing discrimination between groups of items is 
smaller than that within groups of items can be rejected. A small significance level implies a 
very significant rejection of this hypothesis, in turn implying high confidence in the hypothesis 
that the aggregate F-distributed ratios (for the subset) representing discrimination between groups 
of items is greater than or equal to that within groups of items. In this case, the optimal subset 
iS of the component vectors is selected by identifying the subset of component vectors whose 
probability value for the statistic y/ s associated with the subset is a minimum. This minimum 
probability value implies a maximum confidence in the hypothesis that the aggregate 
F-distributed ratios (for the subset) representing discrimination between groups of items is 
greater than or equal to that within groups of items. Selection of the optimal subset S is 
preferably accomplished as follows. 

In step 211, the fk values generated in step 207 are ranked in descending order (i.e., 
from largest to smallest). 

In step 213, a loop is performed over subsets of component vectors whereby, for each 

subset, the statistic ¥ ^associated with the subset is generated and a probability value for statistic 
y/ s is calculated. Preferably, the loop is performed over a counter X ranging from 1 to Z). In 
each iteration of the loop, the following operations are performed. First, the largest X fk values 
identified in step 21 1 are added together as follows 

(Fx = y/ s + fh + —fkx) ) > anc * the resultant sum is normalized by a 
division by X (i.e., F x = F x I X). Second, a probability value, denoted p x , for the normalized 
sum Fx is calculated. Third, the probability value p x and associated counter X are stored. 

The probability value p x for a normalized sum Fx may be calculated in step 213 as 
follows: 

p x = Q(F X | N-1 9 X) 

A more detailed description of the calculation of the probability value p x for a subset of 
component vectors may be found in Press et al., "Numerical Recipes", Cambridge University 
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Press, 1986, pp. 175-190, herein incorporated by reference in its entirety. And a more detailed 
description of the F distribution function is set forth in Freund, "Mathematical Statistics", 5th 
Ed., Prentice Hall, 1992, pp. 314-315, herein incorporated by reference in its entirety. 



In step 215, one or more probability values p x generated in step 213 that satisfy a 
predetermined criterion are selected, and the corresponding subset of component vectors 
corresponding to selected probability value(s) is identified as the optimal subset S of component 
vectors. Preferably, in step 215, the minimum probability value p x generated in step 213 is 
selected, and the corresponding subset of component vectors corresponding to selected 
probability value is identified as the optimal subset S of component vectors. 

Finally, in step 217, one or more of the descriptor vectors for the series of items are 
mapped to a space corresponding to the optimal subset S of component vectors. This mapping 

is preferably accomplished for a given descriptor vector Xy by performing a loop over each 

component vector w belonging to the optimal subset S of component vectors whereby, in each 

iteration of the loop, the contribution w xy (where w is the transpose of the unit vector w for 

the given component vector ~w) is added to a running sum. Consider an example wherein the 

— * -» — » 

optimal subset S includes three component vectors ( w l , w 2 9 w 3 ) . In this example, the 
mapping, denoted M, may be represented by the following: 



FIG. 3 illustrates an example of the operation of the computer processing system 100 in 
determining the set of component vectors (i.e., those w) that maximize the criterion function 



M(X0) = O 



' T 



) W) + (w 2 T Xy)w2 + (w 3 r xfj) W3 
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f(w). The approach defines a vector y that is a function of w along which the contribution from 
the denominator of the criterion function f(w) is independent. 

In step 301, the set of eigenvalue/eigenvector pairs of the matrix e w are calculated. This 
may be accomplished, for example, with the techniques set forth in Press et al., "Numerical 
Recipes", Cambridge University Press, 1986, pp. 349-363, herein incorporated by reference in its 
entirety. The set of eigenvalue/eigenvector pairs for which the eigenvalue is non-zero, denoted 
set E w , is then stored by the computer processing system 1 00. For the sake of description, the 
number of eigenvalue/eigenvector pairs belonging to set E w is denoted K, and a given 
eigenvalue/eigenvector pair belonging to set E w is denoted fa , Vk )■ 

The operations may then define a vector, denoted y , as a function of w as follows: 
y =(ew) 2 w 

where (Sw) 1 =S 4 Vk fe) 1 

where Vk is the unit vector corresponding to eigenvector Vk and 
the 2 operation sums over the K eigenvalue/eigenvector pair 
belonging to set E w . 

The criterion function/(M>) may be rewritten as f(y) as follows: 

f(y)= f teJ~ T ZbiswT 1 y 

where y is the unit vector corresponding to the vector ~y . 

The set of vectors ~y that maximize the criterion function f(/y) may be generated by 
solving the following eigenvalue equation: 
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Finally, the set of vectors y that maximize the criterion function f(y) may then be 
transformed to generate the corresponding set of component vectors w that maximize the 
criterion function/ (w) as follows: 

w = (£ w )"2 y 

In step 303, the computer processing system 100 calculates the matrix 
(fiw) 2 efc(fiw) 2 utilizing the eigenvalue/eigenvector pairs belonging to set E w calculated in 
step 301. Note that (fiw)^ may be calculated as follows: 

(eJ 2 =Z 4 v/fo) 1 

where V/t is the unit vector corresponding to eigenvector vt and 
the 2 operation sums over the K eigenvalue/eigenvector pair 
belonging to set E w . 

In addition, in step 303, the computer processing system 100 calculates the set of 

eigenvalue/eigenvector pairs for the matrix (g w ) 2 £b(£w) 2 . This may be accomplished, for 
example, with the techniques set forth in Press et al., "Numerical Recipes", Cambridge 
University Press, 1986, pp. 349-363, incorporated by reference above in its entirety. The set of 
eigenvalue/eigenvector pairs for which the eigenvalue is non-zero, denoted set E wb , is then stored 
by the computer processing system 100. For the sake of description, the number of 
eigenvalue/eigenvector pairs belonging to set £U is denoted K \ and a given 

eigenvalue/eigenvector pair belonging to set E wb is denoted ( , )• The K 1 eigenvectors 
(each denoted v^) of the eigenvalue/eigenvector pairs belonging to set E w b represent the set of 
vectors ~y that maximize the criterion function fCy^). 
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In step 305, the computer processing system 100 transforms the K ' eigenvectors (each 
denoted v^') of the eigenvalue/eigenvector pairs belonging to set E wb to a set of corresponding 
component vectors in ~w space (each denoted >vt ) as follows: 



The K ' set of vectors (each denoted w k ) represent the set of component vectors that maximize 
the criterion function f(~w\ 



As described above, the computer implemented method of the present invention provides 
a mapping (i.e., a transformation) for the descriptor vectors for a series of items to a space that 
optimally discriminates between groups of items. The method may be used in many domains. 

For example, in the domain of bio-informatics, the items of interest may be genotypes 
that are partitioned into groups based upon phenotypes exhibited by such genotypes; and the 
descriptor vectors associated with such genotypes may represent biological, chemical, and/or 
physical properties of such genotypes. Typically, a candidate genotype and associated descriptor 
vector is identified and one wishes to ascertain to which group the candidate genotype belongs. 
The method discussed above may be used to provide a suggestion as to which group the 
candidate genotype belongs. 

More specifically, the method described above is used to map the descriptor vector for 
each genotype to a space that optimal discriminates between groups. In addition, the statistical 
mean of the mapped descriptor vectors for each group [or some other statistical variable (such as 
the covariance about the mean) based upon the mapped descriptor vectors] is calculated. In 
addition, the mapping function that optimally discriminates between groups is applied to the 
descriptor vector for the candidate genotype. A suggestion as to which group the candidate 
genotype belongs is then determined based upon differences between the mapped descriptor 




— > 
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vector for the candidate genotype and the statistical mean of the mapped descriptor vectors for 
the groups. 

In order to illustrate this example in more detail, consider the case where the genotypes 
are a population of male patients partitioned into two groups (A,B) based upon prostate cancer 
activity experienced by the patients (i.e., patients belonging to group A have not experienced any 
prostate cancer activity, and patients belonging to group B have experienced prostate cancer 
activity). In addition, the descriptor vector associated with each male patient is expression data 
acquired from a gene probe array system, such as the GeneChipsystems developed by 
Affymetrix, Inc. (information is available at http://www.affymetrix.com/products/). In addition, , 
a candidate and associated descriptor vector has been identified. Again, the descriptor vector 
associated with the candidate is expression data acquired from a gene probe array system, and 
one wishes to ascertain to which group the candidate belongs. The method discussed above may 
be used to provide a suggestion as to which group the candidate belongs. More specifically, the 
method described above isVised to map the descriptor vector for each patient to a space that 
optimal discriminates between groups. In addition, the statistical mean of the mapped descriptor 
vectors for each group (or some other statistical variable, such as the covariance about the mean, 
based upon the mapped descriptor vectors) is calculated. In addition, the same mapping function 
is applied to the descriptor vector for the candidate. Finally, a suggestion as to which group the 
candidate belongs is determined based upon differences between the mapped descriptor vector 
for the candidate and the statistical mean of the mapped descriptor vectors for the groups A,B. 

In another example in the domain of bio-informatics, the items of interest may be gene 
sequences that are partitioned into groups; and the descriptor vectors associated with such gene 
sequences may represent expression activities of such gene sequences. The method discussed 
above may be used to derive a probability value for the statistic y/ s that represents a significance 
level at which an hypothesis that the aggregate F-distributed ratios representing discrimination 
between groups of expression activities of gene sequences is smaller than that within groups of 
expression activities of gene sequences can be rejected. A small significance level implies a very 
significant rejection of this hypothesis, in turn implying high confidence in the hypothesis that 
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the aggregate F-distributed ratios representing discrimination between groups of expression 
activities of gene sequences is greater than or equal to that within groups of expression activities 
of gene sequences. In the event that this probability value is less than a threshold value, there is a 
suggestion that the gene sequences express differences in activities. 

In the domain of computational biology and chemistry, for example, the items of interest 
may be molecular complexes (a molecule or portion of a molecule) that are partitioned into 
groups based upon a biological activity (for example, molecular complexes belonging to a group 
all bind to another molecular complex); and the descriptor vectors associated with such 
molecular complexes may represent the physical property (structure, charge distribution), 
biological property and/or chemical property of such molecular complexes. For instance, the 
descriptor vectors associated with a molecule may be based upon data representing the structure 
and/or charge distribution of molecule as outlined in i) U.S. Patent No. 5,784,294 to Piatt et al., 
which is commonly assigned to the assignee of the present invention, ii) R.D. Cramer et al., 
"Comparative Molecular Field Analysis (CoMFA). Effect of Shape on Binding of Steroids to 
Carrier Proteins, J. Am. Chem. Soc. Vol. 110, 1988 pp. 5959-5967, iii) A.C, Good et al, 
"Structure- Activity Relationships from Molecular Similarity Matrices," J. Med. Chem, Vol. 36, 
1993, pp. 433-438, iv) A. Jain et al, "Predicting Biological Activities from Molecular Surface 
Properties. Performance Comparisons on a Steroid Benchmark, " J. Med. Chem, Vol. 37, 1994, 
pp. 2315-2327, v) W. Fisanick et al, "Similarity Searching on CAS Registry Substances, 1: 
Global Molecular Property and Generic Atom Triangle Geometric Searching," Journal of 
Chemical Information and Computer Sciences, Vol. 32, No. 6, 1992, pp. 664-674, and W. 
Fisanick et al, "Similarity Searching on CAS Registry Substances, 2: 2D Structural Similarity," 
Journal of Chemical Information and Computer Sciences, Vol. 34, No. 1, 1994, pp. 130-140, 
hereinafter incorporated by reference in their entirety. Similarly, descriptor vectors may be 
associated with a portion of a molecule (i.e., a subset of the atoms that make Up a molecule) and 
represent the structure and/or charge distribution of the portion of the molecule. In such 
domains, typically, a candidate molecular complex and associated descriptor vector is identified 
and one wishes to ascertain to which group the candidate molecular complex belongs. The 
computer implemented method of the present invention described above may be used to provide 
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a suggestion as to which group the candidate molecular complex belongs. More specifically, the 
method described above is used to map the descriptor vector for each molecular complex to a 
space that optimal discriminates between groups. In addition, the statistical mean of the mapped 
descriptor vectors for each group [or some other statistical variable (such as the covariance about 
the mean) based upon the mapped descriptor vectors] is calculated. In addition, the mapping 
function that optimally discriminates between groups is applied to the descriptor vector for the 
candidate molecular complex. A suggestion as to which group the candidate molecular complex 
belongs is then determined based upon differences between the mapped descriptor vector for the 
candidate molecular complex and the statistical mean of the mapped descriptor vectors for the 
groups. 

The description above illustrates application of the present invention in one-way analysis 
of variance. However, the present invention is not limited in this respect, and can be applied to 
multi-factor designs with multi-way analysis of variance. An example of such a design for an 
item of interest is illustrated in FIG. 4. In this design, a series of items are classified into N 
groups ( N is an integer greater than 1), wherein each group is identified by an identifier i 
ranging from 1 to N represented by the columns of the matrix. A series of factors (typically each 
factor represents one or more experimental treatments) are attributed to the items belonging to 
the groups, wherein each factor is identified by an identifier j ranging from 1 to M represented 

by the rows of the matrix. An identifier Wlyk identifies item k that is attributed to the 
group/factor pair ij , and vector X~ijk identified a descriptor vector attributed to the item 

m ijk . Preferably, the descriptor vectors are stored as part of a file in persistent storage and 
loaded into non-persistent storage for use by the CPU 103 as needed. 

In such a system, one may break any given descriptor vector Xyk into the following 
components: 

— > -» -» — > 

Xyk = a+fii + Jj +Sij +Eijk , where 
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a represents the mean of the descriptor vectors for all the items (of all the groups and all 
the factors); 

fit represents the mean of the descriptor vectors for all the items in group / (of all the 
factors); 

^represents the mean of the descriptor vectors for all the items in factor j (of all the 
groups); 

Sy represents the mean of the deviation of the descriptor vectors for all the items in 

-> ~> -> 

group/factor pair ij from the sum ( a +/?/ + Jj)\ and 

represents deviation from the mean dy to the descriptor vector Xyk . 



These values may be calculated as follows. For each group/factor pair ij of items (i.e., 
/ ranging from 1 to N 9 and j ranging from 1 to M), the mean of the descriptor vectors within 

group/factor pair/, j , denoted X/^ , may be calculated as follows: 

~Xy. - Tig S X~ijk where Yl\j represents the number of items in the 

group/factor pair ij , and the X operator sums over 

the range (k = 1 .. ftij) of items in the group/factor 
pair ij 

Then, for each group / of items (i.e., / ranging from 1 to AO, the mean xT. may be calculated as 
follows: 



x l. - rii 2jXy. where n i represents the number of items in the 
group i , and the 2 operator sums over 
the range (1 M) of factors 

Then, for each factor j (i.e., i ranging from 1 to M) 9 the mean x j. may be calculated as follows: 
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xj.= nj YdXy. where «/ represents the number of items in the 
factor j , and the 2 operator sums over 
the range (1.. iV) of groups 



Then, the mean ~a may be calculated as follows: 

~a= ~n 2 xT where N represents the number of groups , and 
the £ operator sums over the range (1.. N) of 
groups; or equivalently 



~a = ~M^xJ. where M represents the number of factors, and 
the 2 operator sums over the range (1 .. M) of 
factors 



For each group i of items (i.e., i ranging from 1 to N), the mean P> may be calculated as 
follows: 

^> _ -> 

Pi =*/.. - « 



For each factor j of items (i.e., j ranging from 1 to M), the mean ft may be calculated as 



follows: 



For each group/factor pair ij of items (i.e., / ranging from 1 to N, and; ranging from 1 to M), 



the mean (5y may be calculated as follows: 



- */;. - X./ 
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Finally, for each descriptor vector in the group/factor pair ij\ the deviation Syk may be 
calculated as follows: 

The total co variance, denoted s , may be represented by the sum of the following four (4) 

terms: 

(1) 2^ M" fii fit , where the 2 operator sums over the range of groups 

i = 1 to N and M is the number of items in group i. 

(2) 'jj'yf , where the £ operator sums over the range of factors 

and Nj is the number of items in factor i. 

(3) AjNy Si Si , where the 2 operator sums over the range of 

group/factor pairs i,j and Ny is the number if items 
belonging to the group/factor pair i,j. 
— > — > T ^ 

£ijk £ijk , where the 2* operator sums over the range of 



group/factor pairs i j and k items in each group/factor pair 



Note the first term represents the covariance between groups of items; the second term represents 
the covariance between factors; and the third term represents the covariance of the interaction 
between the groups and factors. 

Generally, the principles of the present invention as described above may then be used to 
map the components of the descriptors vectors of the items to a space that: 

1) optimally discriminates between the groups of items; 

2) optimally discriminates between factors; or 

3) optimally discriminates interactions between the groups and factors. 
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More specifically, in order to map the components of the descriptors vectors of the items 
to a space that optimally discriminates between groups of items (case 1 above), an F distributed 
criterion function is used to determine a set of component vectors that maximize the criterion 
function . The criterion function has a numerator and a denominator, whereby the numerator is 
based upon the co variance between groups of items (the first term described above) and the 
denominator is based upon the fourth term described above. The criterion function has the 
general form: 



The set of component vectors is determined by solving for those w that maximize the criterion 
function and an optimal subset of the component vectors is identified. A more detailed 

description of this operation is described above with respect to steps 205- 215 of FIG. 2. Finally, 
the components of the descriptor vectors for the items are mapped to a space corresponding to 
the optimal subset of the component vectors. The resultant data optimally discriminates between 
groups of items. A more detailed description of this operation is described above with respect to 
step 217 of FIG. 2. 

In order to map the components of the descriptors vectors of the items to a space that 
optimally discriminates between factors (case 2 above), an F distributed criterion function is used 
to determine a set of component vectors that maximize the criterion function . The criterion 
function has a numerator and a denominator, whereby the numerator is based upon the 
covariance between factors (the second term described above) and the denominator is based upon 
the fourth term described above. The criterion function has the general form: 




where w is some vector 




-» where w is some vector 
w J 
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The set of component vectors is determined by solving for those w that maximize the criterion 
function /("w), and an optimal subset of the component vectors is identified. A more detailed 
description of this operation is described above with respect to steps 205- 21 5 of FIG. 2. Finally, 
the components of the descriptor vectors for the items are mapped to a space corresponding to 
the optimal subset of the component vectors. The resultant data optimally discriminates between 
factors. A more detailed description of this operation is described above with respect to step 217 
of FIG. 2. 

In order to map the components of the descriptors vectors of the items to a space that 
optimally discriminates interactions between the groups and factors (case 3 above), an F 
distributed criterion function is used to determine a set of component vectors that maximize the 
criterion function . The criterion function has a numerator and a denominator, whereby the 
numerator is based upon the covariance of the interaction between the groups and factors (the 
third term described above) and the denominator is based upon the fourth term described above. 
The criterion function has the general form: 



w 



w 



where w is some vector 



The set of component vectors is determined by solving for those w that maximize the criterion 
function /( >v)» and an optimal subset of the component vectors is identified, A more detailed 
description of this operation is described above with respect to steps 205- 215 of FIG. 2. Finally, 
the components of the descriptor vectors for the items are mapped to a space corresponding to 
the optimal subset of the component vectors. The resultant data optimally discriminates 
interaction between the groups and factors. A more detailed description of this operation is 
described above with respect to step 217 of FIG. 2. 
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The application of the present invention to multi-factor designs may be used in many 
domains. 

In the domain of insurance, for example, the items of interest may be individuals that are 
partitioned into groups based upon certain characteristics such as age, gender, etc. Categories of 
auto insurance polices (for example, categorized based upon the limits/deductable amount for 
liability coverage and/or collision coverage) attributed to the groups of individual may represent 
the factors of the design. The descriptor vectors associated with such individuals may represent 
risk of the individual (for example, may be a dollar amount of the moneys paid to the individual 
arising from auto insurance claims during a predetermined period of time). In this example, the 
computer implemented method discussed above may be used to provide a mapping of the 
components of the descriptor vectors to a space that optimal discriminates between groups of 
individuals; or optimally discriminates between the policies, or optimally discriminates 
interaction between the groups and policies. 

In the domain of agriculture, for example, the items of interest may be plant species that 
are partitioned into groups based upon certain characteristics such as genetic makeup of the plant 
species. Categories of fertilizers (or pesticides) that are applied to the groups of plant species 
may represent the factors of the design. The descriptor vectors associated with such plant species 
may represent a characteristic of the plant species such as yield and/or drought resistance. In this 
example, the computer implemented method discussed above may be used to provide a mapping 
of the components of the descriptor vectors to a space that optimal discriminates between groups 
of plant species; or optimally discriminates between the fertilizers (or pesticides), or optimally 
discriminates interaction between the groups and fertilizers (or pesticides). 

It should be noted that the information encoded by the descriptor vectors may cause 
membership in a group associated with an item of interest. For example, the behavior of a 
molecular complex, which is encoded by a descriptor vector associated with the molecular 
complex, may cause membership in a group (or category) of molecular complexes (for example, 
hydrophilic, polar, "active" with respect to a class of reactions, etc.). In the alternative, the 
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information encoded by the descriptor vector associated with an item of interest may be a 
response to group membership. For example, the risks associated with an insurance policy may 
be encoded by a descriptor vector associated with the insurance policy, whereby the risks are 
dictated by the group (i.e., type) membership of the policy. 

In addition, the present invention may be used to identify molecules within a large 
database of molecules that are similar to a query molecule. The methodology may be 
conceptually divided into three phases: a training phase, an acquisition phase, and a recognition 
phase. 

In the training phase, for each molecule in a series of molecules (denoted the training 
set), the following steps are performed for a given molecule in the training set. Data that 
represents the atomic structure of the given molecule in the training set is stored in persistent 
memory. In addition, data that represents one or more properties of the given molecule is stored 
in persistent memory. The data that represents the atomic structure of the given molecule is read 
from persistent memory and used to define a set of three-dimensional regions of space that 
contain portions of the given molecule. For each region of the given molecule, one or more 
properties of the given molecule are generated and mapped to a one or more property values for 
the grid points of the region; and the properties values for the grid points of the region are used to 
determine a descriptor vector associated with the region that characterizes the region. Preferably, 
one or more components of the descriptor vector associated with a given region represent 
position and/or orientation of axes derived from the property distribution of the region, wherein 
the axes are invariant with respect to translation and rotation of the region. In addition, a group 
identifier is assigned to the region. The group identifier identifies the group to which the region 
belongs. For example, the regions may be partitioned into groups based upon the charge 
distribution of the region (e.g., charged, neutral, polar, non-polar, etc.) or other property or 
behavior of the region. A mapping of one or more components of the descriptor vector for the 
regions is then calculated (preferably by applying the method discussed above is used to generate 
a mapping of components of the descriptor vectors to a space that optimally discriminates 
between the groups of regions - step 217 of FIG. 2). The mapping generated in the training phase 
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is stored in persistent storage and used in the acquisition phase and recognition phase. A more 
detailed description of the training phase for a particular region is illustrated in FIG. 5. 

In the acquisition phase, for each molecule in a series of molecules (which is typically 
distinct from the series of molecules processed in the training phase), the following steps are 
performed for a given molecule. Similar to the training phase, data that represents the atomic 
structure of the given molecule is used to define a set of regions that contain portions of the 
given molecule, and descriptor vectors are generated for such regions of the given molecule. 
Similar to the training phase, one or more items of the descriptor vector associated with a given 
region preferably represents position and/or orientation of axes derived from the property 
distribution of the region, wherein the axes are invariant with respect to translation and rotation 
of the region. The mapping generated in the training phase is used to map one or more 
components of the descriptor vector for each given region (preferably to a space that optimally 
discriminates between groups of regions), and a key is generated based upon the mapping of the 
component(s) of the associated descriptor vector. The key identifies an entry in a multi-map 
(described below) stored in persistent memory. Data identifying the given region, data (or a 
pointer to such data) characterizing a set of axes derived from the property distribution of the 
region, and preferably other data (for example, data identifying the molecule to which the region 
belongs, and data representing the geometric center of the molecule to which the region belongs) 
are then stored in the multi-map at a location identified by the key. As described below in more 
detail, the data characterizing a set of axes derived from the property distribution of a given 
region preferably characterizes transformation between an input reference frame and the inertial 
reference frame for the given region. A more detailed description of an exemplary embodiment 
of the acquisition phase for a particular region is illustrated in FIG. 6. Note that if there are 
common molecules (or regions) in the series of molecules processed in the training phase and 
acquisition phase, data generated in the training phase may be used in the acquisition phase, and 
the operations that use and/or generate such data may be bypassed accordingly. 

In the recognition phase, a query molecule is provided. Similar to the training and 
acquisition phase, data that represents the atomic structure of the query molecule is used to 
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define a set of regions that contain portions of the query molecule, and descriptor vectors are 
generated for such regions of the query molecule. Similar to the training phase, one or more 
items of the descriptor vector associated with a given region preferably represents position and/or 
orientation of axes derived from the property distribution of the region, wherein the axes are 
invariant with respect to translation and rotation of the region. For each query region, the 
mapping generated in the training phase is used to map one or components of the descriptor 
vector for the given query region (preferably to a space that optimally discriminates between 
groups of regions), and a key is generated based upon the mapping of the component(s) of the 
associated descriptor vector. For each key, a multi-map entry identified by the key is retrieved 
and data stored therein [i.e., data identifying one or more regions associated therewith, data 
characterizing a set of axes derived from the property distribution of each region associated 
therewith, data identifying the molecule to which each region associated therewith belongs, and 
data representing the geometric center of the molecule to which each region associated therewith 
belongs] are read from the table. For each region identified by the retrieved table entry, an 
hypothesized match is constructed and added to a vote table. After processing all of the regions 
identified by the retrieved table entry, the vote table is sorted to determine a set of potential 
matching regions, and the set of potential matching regions (and/or the molecules to which the 
set of potential matching regions belong) is made accessible to the user via an I/O device. A 
more detailed description of an exemplary embodiment of the recognition phase for a region of 
the query molecule is illustrated in FIGS. 7(A) and (B). 

Training Phase 

In the training phase, for each molecule in a training set, data that represents the structure 
of the given molecule is stored in persistent memory. Preferably, the data represents the atomic 
structure of the given molecule in an arbitrary three-dimensional reference frame, which is 
referred to below as the "input reference frame". The data may be obtained from a database 
(public or private) or be derived by traditional molecular modeling techniques. 
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In addition, for each molecule in the training set, data that represents the property 
distribution of one or more properties of the given molecule is generated, for example by reading 
such data from persistent memory. The property distribution for the property of a given molecule 
may be an atomic property distribution (i.e., data representing a heuristic, such as 
electronegativity, hydrophobicity or polarity, that characterizes the behavior or property of an 
atom in the given molecule), a surface property distribution (i.e., data representing a heuristic that 
characterizes the behavior or property of a surface on or within the molecule), or a volumetric 
property distribution (i.e., data characterizing the behavior or property of a region of the 
molecule- for example, electron density). 

In addition, for each given molecule in the training set, the structure data is used to 
generate a set of three-dimensional regions of space, hereinafter referred to "scoops' 1 , in the input 
reference frame. Preferably, each scoop is a spherical region in the input reference frame having 
a center (i.e., a point in the input reference frame) and a radius. Preferably, the center of each 
scoop corresponds to one or more heuristics of the given molecule. For example, the heuristic 
may correspond to the coordinate of i) the nucleus of one or more atoms of the given molecule, 
ii) the center of one or more bonds between atoms of the given molecule, iii) the end point of an 
extension to a bond (typically the length of the extension is a factor of bond length) between 
atoms of the given molecule, iv) a grid points in the given molecule, v) the geometric center of 
the molecule, vi) ring centers in the given molecule, or vii) lone pairs of electrons in the given 
molecule. In addition, the radius of each scoop may be set to a predetermined value (for 
example, 3 angstroms, or 5 angstroms) or may be based upon heuristics of the given molecule. 
It should be noted that a scoop can be any arbitrary three-dimensional region in the input 
reference frame. 

As illustrated in FIG. 5, the training phase preferably performs a nested loop over the 
molecules in the training set; wherein, for a given molecule, each scoop in the set of scoops for 
the given molecule is processed as follows. 
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In step 501, the property distribution of one or more properties of the given molecule that 
is relevant to the scoop are identified, and the relevant property distribution is mapped to a 
property field that represents the value of the property at points in the input reference frame. 
There are many well-known techniques to map a property distribution (atomic, surface, 
volumetric) to a property field, the specifics of which is not relevant to the present invention. 
Examples of such techniques may be found in U.S. Patent 5,025,388 to Cramer, III et al., 
hereinafter incorporated by reference in its entirety. With regard to atomic properties for the 
given molecule, preferably the properties of atoms that are contained within the scoop are 
identified as relevant properties for the scoop in step 501. With regard to surface properties for 
the given molecule, preferably the properties of surfaces that are contained (partially or fully) 
within the scoop are identified as relevant properties for the scoop in step 501 . With regard to 
volumetric properties for the given molecule, preferably the properties of volumes that are 
contained (partially or fully) within the scoop are identified as relevant properties for the scoop in 
step 501. 

In step 503, the property field generated in step 501 is mapped to grid points contained in 
the scoop to determine contribution of property field at the grid points. Preferably, the grid 
points are evenly spaced within the given scoop. In addition, the contribution of the property 
field at a given grid point preferably includes two values: the first value is a positive value, 
which for the sake of description is referred to below as a mass value; the second value is a real 
value, which for the sake of description is referred to below as a charge value. There are many 
possible techniques to calculate the contribution of a property field at a given grid point, the 
specifics of which is not relevant to the present invention. For example, the following "smearing 
function" can be used in step 503 to map values of the points of the property field to a first value 
jUi (mass value) at grid point /: 




where Aj represents the value of the property field at a point j in the 
property field, 

Ry represents the distance between the point j in the property field 
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and the point / in the scoop; 
Lo is a damping factor, which typically ranges between .5 and 
2.5; 

B is an arbitrary constant, which typically is 1 ; and 

2 sums over all the points in the property field for the scoop. 

Alternate smearing functions may be gaussian forms. In addition, the following equation may be 
used to map values of the points of the property field to a second value Xi (charge value) at grid 
point i: 

Xi 

where Ji represents the mean first value (mass value) for all the grid 
points of the scoop 

In step 505, a descriptor vector for the scoop is generated based upon the property values 
of the grid points of the scoop. Preferably, the components of the descriptor vector for the scoop 
includes one or more of the following data values, details of which are set forth below: 

ID data value identifying the molecule from which the scoop is 

derived 

N an integer value representing the number of grid points for the 

scoop 

M total mass for the scoop, which is a sum of the first values (i.e., 

"mass values" or fit) for the grid points of the scoop 

Q total charge for the scoop, which is the sum of the second values 

(i.e., "charge values" or Xi ) f° r the grid points of the scoop 
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Sx ? Sy , S Z 



components of a vector describing translation between origin of 
the input reference frame and center of mass in the input reference 
frame 



I x 9 I y , I z principal values of the moment of inertia tensor of the scoop with 

respect to the center of mass of the scoop 

Se, So> , Sn So> , Se are polar angles that describe an axis of rotation. 

Sn represents an angle of rotation about that angle. 
Together, the angles represent a rotation transformation between 
principal axes of the moment of inertia tensor and axes of the input 
reference frame 



v x , v y , v z components of a vector describing translation between the center 

of scoop in the inertial reference frame and center of mass in the 
inertial reference frame 



d x , d y , d z components of a vector describing translation between center of 

mass in the inertial reference frame and center of dipole in the 
inertial reference frame 



Ci x , Ci y , Ci z components of a vector describing third order moment of mass 
about a center of expansion in the internal reference frame 

Cq x , Cq y , Cq z components of a vector describing third order moment of charge 
about a center of expansion in the internal reference frame 

Qxx , Qxy , Qxz components of a tensor characterizing quadrupolar moment about a 
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Qyy > Qyz 9 Qzz center of expansion in the internal reference frame 

Note that the components ID, N, M 5 Q, I x , I y , I z , v x , v y , v z , d x , d y , d 2 , Ci x , Ci y , Ci z , 
Cq x , Cq y , Cq z , Q xx , Q xy , Q xz , Qyy , Q yz , Qzz characterize the property distribution with 
respect to invariant axes (i.e., axes that are invariant with respect to translation and rotation of the 
given scoop - for example the sensed inertial axes of the given scoop). Also note that the 
components S x , S y , S z , Se 9 S<p , Sn are variant with respect to translation and rotation of the 
scoop, and characterize orientation of the invariant axes with respect to the input reference frame 
of the scoop. A more detailed description of the calculation of the components of the descriptor 
vector for a given scoop is set forth below. 

In step 507, a group identifier is assigned to the scoop. The group identifier identifies the 
group to which the region belongs. For example, the scoops may be partitioned into groups 
based upon the charge distribution of the scoop (e.g., charged, neutral, polar, non-polar, etc.) or 
other property or behavior of the scoops. 

In step 509, statistical analysis is used to generate a mapping of one or more components 
of the descriptor vector for the scoop. Preferably, the computer implemented method of the 
present invention discussed above is used to generate a mapping of the components of the 
descriptor vector for the scoop that characterize position and/or orientation of invariant axes 
derived from the property distribution of the given scoop (i.e., ID, N, M, Q, I x , I y , I z , v x , V y , 
v z , d x , d y 5 d z , Ci x , Ci y , Ci z , Cq x , Cq y , Cq z , Q xx , Q xy , Q xz , Qyy , Q yz , Qzz) to a 
space that optimally discriminates between the groups of scoops (described above with respect to 

step 217 of FIG. 2). As described above with respect to step 217 of FIG. 2, the mapping is based 

— > — > — > 

upon the transpose of unit vectors for the component vectors w i , w 2 , w 3 . The mapping 
generated in step 509 is preferably stored in persistent storage for subsequent use in the 
acquisition phase and recognition phase. 
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Acquisition Phase 



In the acquisition phase, for each molecule in a series of molecules, data that represents 
the structure of the given molecule is stored in persistent memory. Preferably, the data 
represents the atomic structure of the given molecule in an arbitrary three-dimensional reference 
frame, which is referred to below as the "input reference frame". The data may be obtained from 
a database (public or private) or be derived by traditional molecular modeling techniques. 

In addition, for each molecule in a series of molecules, data that represents the property 
distribution of one or more properties of the given molecule is generated, for example, by 
reading such data from persistent memory. The property distribution for a property of the given 
molecule may be an atomic property distribution, a surface property distribution, or a volumetric 
property distribution as described above. 

In addition, for each given molecule in the series of molecules, the structure data is used 
to define a set of scoops in the input reference frame. Preferably, each scoop is a spherical region 
in the input reference frame having a center and radius as described above. 

As illustrated in FIG. 6, the acquisition phase preferably performs a nested loop over the 
series of molecules, wherein, for a given molecule, each scoop in the set of scoops for the given 
molecule is processed as follows. 

In step 601, properties of the given molecule that are relevant to the scoop are identified, 
and the relevant properties are mapped to a property field that represents the value of the property 
at points in the input reference frame. This operation is similar to the processing described above 
with respect to step 501 of the training phase. 

In step 603, the property field generated in step 601 is mapped to grid points contained in 
the scoop to determine contribution of property field at the grid points. Preferably, the 
contribution of property field at a given grid point includes two values: a positive mass value 
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and a real charge value. This operation is similar to the processing described above with respect 
to step 503 of the training phase. 

In step 605, a descriptor vector for the scoop is generated based upon the property values 
of the gird points of the scoop generated in step 603. Preferably, the components of the 
descriptor vector for the scoop includes one or more of the data values as described above with 
respect to step 505 of the training phase. 

In step 607, the mapping generated in step 509 of the training phase is used to map one or 
more components of the descriptor vector for the scoop (preferably, the component(s) of the 
descriptor vector is mapped to a space that optimally discriminates between the groups of scoops 
as described above with respect to step 217 of FIG. 2). 

In step 609, a key is generated based upon the mapping of the component(s) of the 
descriptor vector for the scoop in step 607. The key identifies an entry in a multi-map stored in 
persistent memory. The multi-map is an associative memory which permits more than one entry 
stored in the memory to be associated with the same key. A detailed description of a multi-map 
is set forth in D.R. Musser and Atul Saini, STL Tutorial and Reference Guide: C++ 
Programming with the Standard Template Library (Addison- Wesley 1996), herein incorporated 
by reference in its entirety. Preferably, the multi-map is formed from a hash table. A more 
detailed description of a hash table may be found in R. Sedgewick, Algorithms in C++ 
(Addison- Wesley 1992), herein incorporated by reference in its entirety. In the alternative, the 
multi-map container may be formed from a linked list data structure, or a tree structure such as 
an AVL-tree or B* tree as described in D.R. Musser and Atul Saini, herein incorporated by 
reference in its entirety. One skilled in the art will realize that there are many possible 
underlying implementations for the multi-map data structure. 

Finally, in step 611, data identifying the given scoop, data (or a pointer to such data) 
characterizing a set of axes derived from the property distribution of the scoop, and preferably 
other data (for example, data identifying the molecule to which the scoop belongs, and data 



Y0998-112 



-34- 



representing the geometric center in the input reference frame of the molecule to which the scoop 
belongs) are then stored in the multi-map at a location identified by the key generated in step 
609. Preferably, the data characterizing a set of axes derived from the property distribution of a 
given scoop characterizes transformation between an input reference frame and the inertial 
reference frame for the given scoop - the components S x , S y , S z , Se , S<d , Sn of the descriptor 
vector for the scoop as discussed above with respect to step 605. 

An exemplary multi-map entry is illustrated in FIGS. 7(A) and (B). The entry 701 
includes a series of segments 703-1 and 703-2 (two shown) coupled via a link-list data structure. 
Each segment includes an first ID field 711 that stores an identifier for a scoop, a pointer 713 to 
data (or the data itself) characterizing a set of axes derived from the property distribution of the 
scoop, a second ID field 715 identifying the molecule to which the scoop belongs, a pointer 717 
to data (or the data itself) representing the geometric center of the molecule to which the scoop 
belongs, and a pointer 719 to the next segment in the table entry. 

The operations of step 601-61 1 are applied for each scoop in a given molecule. At the 
end of the acquisition phase, the multi-map stores entries each corresponding to one or more 
scoops for the series of molecules being studied. 

Recognition Phase 

FIGS. 8(A) and (B) illustrates the recognition phase for a query molecule. In step 801, 
data that represents the structure of the query molecule in the input reference frame is stored in 
persistent memory. In step 803, data that represents the property distribution of one or more 
properties of the query molecule is generated, for example by reading such data from persistent 
memory. The property distribution of the property for a given molecule may be an atomic 
property distribution, a surface property distribution, or a volumetric property distribution as 
described above. In step 805, the structure data of the query molecule is used to define a set of 
scoops (in the input reference frame) in the query molecule. In steps 807-825, a loop is 



Y0998-112 



-35- 



performed over the set of scoops in the query molecule wherein each scoop is processed as 
follows. 

In step 809, the properties of the query molecule that are relevant to the scoop are 
identified, and the relevant properties are mapped to a property field that represents the value of 
the property at points in the input reference frame of the query molecule. This operation is 
similar to the processing described above with respect to step 501 of the training phase. 

In step 811, the property field generated in step 809 is mapped to grid points contained in 
the scoop to determine contribution of property field at the grid points. Preferably, the 
contribution of property field at a given grid point includes two property values as described 
above: a positive mass value and a real charge value. This operation is similar to the processing 
described above with respect to step 503 of the training phase. 

In step 813, a descriptor vector for the scoop is generated based upon the property values 
(mass values and charge values) of the grid points of the scoop generated in step 811. 
Preferably, the components of the descriptor vector for the scoop includes one or more of the 
data values as described above with respect to step 505 of the training phase. 

In step 815, the mapping generated in step 509 of the training phase is used to map one 
or more components of the descriptor vector for the scoop generated in step 813 (preferably, the 
component(s) of the descriptor vector are mapped to a space that optimally discriminates 
between the groups of scoops as described above with respect to step 217 of FIG. 2). This 
operation is similar to the processing described above with respect to step 607 of the acquisition 
phase. 

In step 817, a key is generated based upon the mapping of the component(s) of the 
descriptor vector generated in step 815. The key identifies an entry in the multi-map stored in 
persistent memory. This operation is similar to the processing described above with respect to 
step 609 of the acquisition phase. 
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In step 819, the multi-map entry identified by the key generated in step 817 is retrieved 
and the data stored therein are read from the retrieved multi-map entry. 

In step 821, for each scoop identified by the retrieved multi-map entry (i.e., each entry 
segment in FIG. 7), an hypothesized match is constructed to determine a set of transformation 
parameters whereby a set of axes derived from property distribution of the query scoop is aligned 
with a set of axes derived from property distribution for each molecular scoop identified by the 
retrieved multi-map entry, which for the sake of description is referred to below as the stored 
scoop. 

In step 823, a label corresponding to transformation parameters of the hypothesized 
match is generated, and the hypothesized match is added to a vote table. The vote table is an 
associative memory of entries keyed by the label. Each entry of the vote table stores i) an 
accumulated score of the number of stored scoops whose hypothesized match corresponds to the 
label, and ii) data identifying those stored scoops whose hypothesized match corresponds to the 
label. Preferably, the vote table is implemented as a map. A detailed description of a map is set 
forth in D.R. Musser and Atul Saini, STL Tutorial and Reference Guide: C++ Programming with 
the Standard Template Library (Addison- Wesley 1996), incorporated by reference above in its 
entirety. The map may be formed from a hash table. A more detailed description of a hash table 
may be found in R. Sedgewick, Algorithms in C++ (Addison- Wesley 1992), incorporated by 
reference above in its entirety. In the alternative, the map may be formed from a linked list data 
structure, or a tree structure such as an AVL-tree or B* tree as described in D.R. Musser and 
Atul Saini. One skilled in the art will realize that there are many possible underlying 
implementations for the map data structure. 

In step 827, after processing all of the stored scoops identified by the retrieved multi-map 
entry (step 819), the vote table is sorted to determine a set of potential matching scoops. For 
example, one or more entries with the highest score may be selected, and those scoops identified 
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by the selected entries of the vote table may be selected as the set of potential matching scoops 
for the query molecule. 

Finally, in step 829, the set of potential matching scoops (and/or the molecules to which 
the set of potential matching scoops belong) is made accessible to the user, for example, via 
graphical user interface 907 of the computer processing system. 

In the preferred embodiment of the present invention, the operations of the recognition 
phase in constructing the hypothesized match, adding the hypothesized match to the vote table, 
and sorting the vote table to report results to the user are illustrated in the flow chart of FIGS. 
9(A) and (B). As described above, each entry of the multi-map generated in step 61 1 of the 
acquisition phase preferably stores: i) data characterizing transformation between the input 
reference frame and the sensed inertial reference frame for a given stored scoop - which may 
include the components S x , S y , S z , S e , S<n , Sn of the descriptor vector for the stored scoop; ii) 
data identifying the molecule to which the stored scoop belongs, denoted the stored molecule; 
and iii) data representing the geometric center in the input reference frame of the stored 
molecule. Similar operations are performed in the recognition phase in constructing the 
hypothesized match as illustrated in the flow chart of FIGS. 9(A) and (B). 

More specifically, in step 901, data characterizing transformation between the input 
reference frame and the sensed inertial reference frame for the query scoop - which may include 
the components S x , S y , S z , Se, S<d , So of the descriptor vector for the query scoop is 
generated. 

In step 903, data characterizing transformation between the input reference frame for the 
stored molecule (denoted "stored frame") and the sensed inertial reference frame of the stored 

AAA 

scoop (U i , U2 , Ui of the stored scoop as described below) and data characterizing 
transformation between the input reference frame of the query molecule (denoted "query frame") 

AAA 

and the sensed inertial reference frame of the query scoop (U \ , U2 , U3 of the query scoop as 
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described below) are used to calculate a transformation (a translation and rotation 
transformation) that aligns the two sensed inertial frames, which also represents a transformation 
from the stored frame to the query frame. Details of an exemplary technique for calculating this 
transformation may be found in K. Turkowski, "Graphics Gems," Academic Press, edited by A. 
Glassner, pgs. 522-532, herein incorporated by reference in its entirety. 

In step 905, the alignment transformation generated in step 903 is applied to the data 
representing the geometric center of the stored molecule, which is preferably retrieved from data 
included in the matching entry of the multi-map, to generate data representing center of the 
stored molecule in the query frame. 

In step 907, it is determined if the center of the stored molecule in the query frame lies 
within the volume of the query molecule. The volume of the query molecule may be calculated 
utilizing various well-known techniques, such as the techniques described in Connelly, M.L., 
"Computation of Molecular Volume" JACS, Vol. 107, 1985, pg. 1118-1124, hereinafter 
incorporated by reference in its entirety. 

If the test of step 907 fails (the center of the stored molecule in the query frame lies 
outside the volume of the query molecule), in step 909 the construction of the match hypothesis 
ends without adding an entry to the vote table, and processing continues to step 909 

However, if the comparison of step 907 is successful, in step 91 1 data representing the 
rotation of the alignment transformation generated in step 903 and data representing center of the 
stored molecule in the query frame are preferably quantized to form a integer pair, and a label is 
generated based upon such data and the identifier of the stored molecule. 

In step 913, if the vote table does not include an entry corresponding to the label 
generated in step 91 1, a new entry associated with the corresponding label is added to the vote 
table. The new entry includes a score field with an initial value (for example, 1), data 
identifying the given stored scoop, and data identifying the alignment transformation generated in 
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step 903. Otherwise (the vote table does not include an entry corresponding to the label 
generated in step 911), the score field of the corresponding entry is incremented and possibly 
data identifying the stored scoop is added to the entry. In addition, the alignment transformation 
data may be updated, for example, to represent the cumulative average alignment transformation. 
The operation then continues to step 909. 

In step 909, one or more entries of the vote table is selected, and the alignment 
transformation corresponding to the selected entry (ies) is applied to the corresponding stored 
molecule (which is identified by the label associated with the vote table entry) to generate an 
alignment of the stored molecule in the query frame. 

Finally, in step 915, the alignment of the stored molecule in the query frame is reported to 
the user via an I/O device and operation of the recognition phase ends. 

It should be noted that the operations described above with respect to FIG. 9 represent a 
preferred embodiment of the present invention. One skilled in the art will realize that the 
operations similar to those described above can be used to construct an hypothesized match 
between a stored scoop and a query scoop based upon any data characterizing a set of axes 
derived from the property distribution of the stored scoop and query scoop, respectively. 

As described above, the training phase is used to define the association criteria between 
query scoops and stored scoops, and keys and the corresponding multi-map data structure capture 
the associations. In an alternate embodiment, any arbitrary selection scheme (for example, one 
based on a distance metric) can be used to associate a query scoop with a stored scoop, which 
leads to construction of a match hypothesis between the query scoop and the associated stored 
scoop. 

The computer processing system 1 00 that implements the present invention may be 
distributed in nature as shown in FIG. 1(B). More specifically, a distributed computer processing 
system comprises more than one CPU 103 (three shown 103-1, 103-2,103-3) with each of these 
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CPUs communicating with one another via message passing utility 113. The message passing 
utility 113 may be implemented via shared memory, a network connection, a high speed switch 
or some other method that allows data to be passed from CPU to CPU. A distributed computer 
processing system is preferably used for the recognition phase of the present invention because of 
the inherently parallel nature of the algorithm. More specifically, the multi-map data structure 
generated in the acquisition mode as described below is preferably partitioned amongst the CPUs 
of the distributed system. The multi-map data structure may also be partitioned amongst the 
various nonvolatile storage devices 108 associated with a given CPU 103. For example, the 
multi-map may be partitioned into nine portions MM11, MM12, MM13, MM21, MM22, MM23, 
MM31, MM32, MM33 among three CPUs 103-1,103-2 and 103-3 and their associated nine 
storage devices 108-11, 108-12 ,108-13 , 108-21,1 08-22, 108-23, 108-31, 108-32, 108-33 as 
shown. In addition, the vote table data structure may be similarly partitioned amongst the CPUs 
of the distributed system. When such a system is used in the recognition phase described above, 
as table entries are retrieved from the multi-map, such table entries are routed via the message 
passing utility 1 13 to the appropriate CPU for accumulation in the proper segment of the vote 
table. In the end, a distributed merge sort is preferably used to collate all of the resulting 
hypothesized matches on a single CPU. 

Derivation of Components of Descriptor Vectors for a Scoop 

As described above, the components of the descriptor vector characterizing a given scoop 
preferably includes S x , S y , S z , which represent the components of a vector describing translation 
between the origin of input reference frame and center of mass in the input reference frame. 
These components may be calculated utilizing well known techniques, including the operations 
described in column 6 of U.S. Patent No. 5,784,294 to Piatt et al., incorporated by reference 
above in its entirety. 

In addition, the components of the descriptor vector characterizing a given scoop 
preferably includes the principal values of the moment of inertia tensor I for the scoop, which 
are denoted I x , I y , I 2 . The inertial reference frame for the scoop is preferably characterized by 
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an origin at the center of mass in the input reference frame, and a set of three axes (denoted U \ , 

A A AAA 

U2 , U3) in the input reference frame. The axes U\, U2 , U3 are unit vectors that point in the 
same direction as corresponding eigenvectors (principal axes) of the diagonalized moment of 
inertia tensor I of the scoop. The principal values of the moment of inertia tensor I of the scoop 
(I x , I y , I z ) are represented by the eigenvalues of the diagonalized moment of inertia tensor I of the 
scoop. The eigenvectors (principal axes) and the corresponding eigenvalues of the diagonalized 
moment of inertia tensor I for the scoop may be calculated utilizing the operations described in 
columns 6-8 of U.S. Patent No. 5,784,294 to Piatt et al., incorporated by reference above in its 
entirety, 

AAA 

Importantly, the three axes U\ , U2 , U3 do not sufficiently denote an internal frame of 
reference for coordinate transformation because the signs of the corresponding eigenvectors are 
ambiguous. The signs must be determined from other information. 

A A 

Accordingly, the present invention provides a procedure for sensing the axes U\ , U2 , 

A 

U3 with a third order moment vector, which is denoted the asymmetric vector, in order to 
denote an internal frame of reference sufficient for coordinate transformation. The asymmetric 

vector has a general form C = (C x , C y , C z ). This general form may characterize a third order 

moment of mass about a center of expansion, denoted above as C, = (Ci x , Ci y , Ci z ). In the 
alternative, the general form may characterize a third order moment of charge about center of 

expansion, denoted above as C q = (Cq x , Cq y , Cq z ). 

In the preferred embodiment of the present invention, the components of the asymmetric 
-> 

vector C = (C x , C y , C z ) are derived as follows: 
C = J p(~x^)\!c\ 2 ~x d 3 x 

where 1c represents a vector between the center of expansion and any arbitrary 
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point in the input reference frame, and 



p(lt^) represents the relevant distribution (i.e., represents the 

mass distribution for the asymmetric vector 

characterizing third order moment of mass, or represents the 

charge distribution for asymmetric vector characterizing 

third order moment of charge) of the scoop at points in the input reference 

frame corresponding to the vector it . 

In systems where the distribution p(l?) of the scoop has discrete values over points in the input 

-> 

reference frame, the components of the asymmetric vector C = (C x , C y , C z ) may be derived as 
follows: 



^ - |->l 2 -> 



XV represents 

reference frame and a grid point / in the input reference frame, 



where xv represents a vector between the center of expansion in the input 



Pi is the mass/charge property at the grid point i (when generating the 
asymmetric vector characterizing third order moment of mass,/?/ is the 
mass value ju ; for the grid point i as described above 
with respect to step 503; and when generating the asymmetric vector 
characterizing third order moment of charge,/?; is the charge value Xi 
for the grid point / as described above with respect to 
step 503); and 

£ sums over the grid points of the scoop. 



With respect to the asymmetric vector C, characterizing third order moment of mass for the 
scoop, preferably the center of expansion is the center of mass of the scoop. With respect to the 
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asymmetric vector Cq characterizing third order moment of charge for the scoop, the center of 
expansion may be the center of charge, center of dipole, or center of quadrupole for the scoop. 

As described above, the components of asymmetric vector C (such as the asymmetric 

-> . — > 

vector C, characterizing third order moment of mass for the scoop, or the asymmetric vector Cq 

A A 

characterizing third order moment of charge for the scoop) are used to sense the axes U \ 5 U2 , 

A 

U3 in order to denote an internal frame of reference sufficient for coordinate transformation. 

AAA 

Preferably, the axes U\ , U2 , U3 are sensed by looping over a counter y ranging from 1 to 2 
with an increment of 1. In each iteration of the loop, the following conditional operation is 

A 

performed: if the dot product of the axis U\y] and the asymmetric vector is less than zero, then 

A A A 

the sign of U\y\ is swapped (i.e., U\y\ = - W[y]). In the first iteration of the loop, the counter y 

A 

is 1 and the conditional operation tests whether the dot product of the axis U \ and the 

A A 

asymmetric vector is less than zero. If so, the sign of U \ is swapped. Otherwise, the sign of U \ 
remains unchanged. In the second iteration of the loop, the counter^ is 2 and the conditional 

A 

operation tests whether the dot product of the axis U2 and the asymmetric is less than zero. If 

A A 

so, the sign of U2 is swapped. Otherwise, the sign of U2 remains unchanged. The loop then 

A A A 

terminates, and the third axis U3 is constructed as the cross product of U\ and W2. 

AAA 

Importantly, the sensed axes U\ , U2 , U3 denote a reference frame suitable for 
coordinate transformation. More specifically, this reference frame is characterized by an origin 

AAA 

at the center of mass in the input reference frame, and the sensed axes ( U \ , U2 , U3 in the input 
reference frame). This reference frame is preferably used to generate components of the 
descriptor vector that characterizes a scoop. 

In addition, the components of the descriptor vector characterizing a given scoop includes 
the components of a transformation matrix representing rotation between the principal axes of 
moment of inertia tensor I for the scoop and the axes of input reference frame. Preferably, this 
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transformation matrix, which for the sake of description is denoted R , is derived from the sensed 
axes u\ (= u\ x x + u\ y y + u\ z z) , ui (= u^x + u 2y y + u 2z z), and w 3 (= ui x x + ui y y + m z z) as 
follows: 



Note that this representation has nine terms. Preferably, the transformation matrix R may be 
represented by three angles: Se, So , Sn. So , Se are polar angles that describe an axis of 
rotation. So represents an angle of rotation about that angle. Together, the angles represent the 
rotation transformation represented by the transformation matrix R. The rotation matrix R and 
the angles Se, So , Sn can be calculated using the technique described in M. Pique, "Graphics 
Gems," Academic Press, edited by A. Glassner, pgs. 465-467, herein incorporated by reference in 
its entirety. Note that this representation of the rotation transformation using the angles Se, So , 
Sn only has three terms, which is advantageous because it lowers the storage allocation 
requirements for the components of the descriptor vector of the scoop that represent the 
transformation. 

In addition, the components of the descriptor vector characterizing a given scoop 
preferably includes the components of a vector describing translation between center of scoop in 
the inertial reference frame and center of mass in the inertial reference frame, which are denoted 
v x , v y , v z . Preferably, the components v x , v y , and v z are derived by generating a vector 
representing translation between center of scoop in the input reference frame and center of mass 
in the input reference frame, and then applying the rotation transformation matrix R described 
above to transform the vector from the input reference frame to the inertial reference frame 

AAA 

represented by the sensed axes U\ 9 U2 , W3. 



U\ X U\y 

R = U2x Uly 
U3x U3y 



U2z 



U3z 
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In addition, the components of the descriptor vector characterizing a given scoop 
preferably includes the components of a vector describing translation between center of mass in 
the inertial reference frame and center of dipole in the inertial reference frame, which are denoted 
d x , d y , d z . Preferably, the components d x , d y , and d 2 are derived by generating a vector 
representing translation between center of mass in the input reference frame and center of dipole 
mass in the input reference frame, and then applying the rotation transformation matrix R 
described above to transform the vector from the input reference frame to the inertial reference 

AAA 

frame represented by the sensed axes U \ , Ui , U3 . The center of dipole in the inertial reference 
frame may calculated utilizing the operations described in columns 10-1 1 of U.S. Patent No. 
5,784,294 to Piatt et. al., incorporated by reference above in its entirety. 

In addition, the components of the descriptor vector characterizing a given scoop 
preferably includes one or more components of a tensor characterizing the quadrupolar moment 
about a center of expansion, which are denoted Q xx , Q xy , Q XZ) Qyy , Q y z , Qzz. Preferably, the 
components Q xx , Q xy , Qxz, Qyy , Q yz , Qzz characterize contribution of such quadrupolar moment 

AAA 

along the sensed axes U \ , U2 , U3 . In this case, components of a tensor Q characterizing 
quadrupolar moment about a center of expansion are generated. The components of the tensor Q 
may be calculated utilizing the operations described in columns 10-11 of U.S. Patent No. 
5,784,294 to Piatt et. al., incorporated by reference above in its entirety. Finally, the following 
operations are performed with respect to the tensor Q to calculate contribution of the tensor Q 

AAA 

along the sensed axes U \ , U2 , U3 : 
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• ♦ 



While the invention has been described in connection with specific embodiments, it will 
be understood that those with skill in the art may develop variations of the disclosed 
embodiments without departing from the spirit and scope of the following claims. 
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